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Galaxy surveys provide a large-scale view of the universe that typically has a lim- 
ited line-of-sight or redshift resolution. The lack of radial accuracy in these surveys 
can be modelled by picturing the universe as a set of concentric radial shells of finite 
width around the observer, i.e, an onion-like structure. We present a new N-body 
simulation with 2048"^ particles developed at the Marenostrum supercomputer with 
the GADGET-2 code. Using the lightcone output we build a set of angular maps 
that mimic this onion-like representation of the universe. The onion maps are a highly 
compressed version of the raw data (i.e., a factor > 1000 smaller size for arcminute res- 
olution maps) and they provide a new and powerful tool to exploit large scale structure 
observations. We introduce two basic applications of these maps that are especially 
useful for constraning dark energy properties: the baryon acoustic oscillations (BAO) 
in the galaxy power spectrum and all-sky maps of the weak lensing distortion. In par- 
ticular, from the matter density maps, we determine the smallest scale where linear 
theory and the Gaussianity of the error analysis applies. Using the weak lensing maps, 
we measure the convergence power spectra and compare it to halo fit predictions. We 
also discuss mass resolution effects and error determinations. As a further application, 
we compute the variance and higher-order moments of the maps. We show that sam- 
pling variance on scales of few degrees is quite large, resulting in a significant (25% at 
10 arminute scales) bias in the variance. We caution that current lensing surveys such 
as the COSMOS HST should take into account this bias and extra sampling error in 
their clustering analyses and inferred cosmological parameter constraints. Finally, we 
test the importance of projection effects in the weak lensing mass reconstruction. On 
the mean, the mass calibration works well but it exhibits a large non-Gaussian scatter 
what could induce a large bias in the recovered mass function. 



1 INTRODUCTION 



Upcoming astronomical surveys will gather many Terabytes 
of unprecedented high quality data containing the relevant 
information to answer key cosmological questions, ranging 
from the nature of the initial conditions in the structure 
formation of the universe, how galaxies and clusters form 
and evolve, or the properties of the so called dark energy 
and theory of gravity on cosmological scales. 

These datasets will pose a great challenge to the sci- 
entific community to develop the appropriate data analysis 
tools to compress the overwhelming raw data into a few 
numbers, eg a set of cosmological parameters. Simulating 
surveys, with their anticipated volume, resolution and com- 
plexity, has become a standard tool to prepare the scien- 



tific exploitation and to understand real astronomical data. 
Thus, analyzing mock surveys suffers from similar limita- 
tions, with the agravant that it requires a large number of 
simulations to pin down statistical errors and explore cos- 
mological parameter space. 

Measuring redshifts for many galaxies is very costly (es- 
pecially at z > 1), even for very large telescopes. Thus, to 
explore the largest scales with catalogs containing many mil- 
lions of galaxies typically requires a photometric approach 
to obtain galaxy redshifts, as is the case for most of the 
upcoming or planned surveys (such as DES, PAU, VHS, 
PanSTARRS, LSST, DUNE). In these surveys information 
in the radial direction is washed out on scales smaller than 
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the photometric error width. This limits the amount of in- 
formation that can be used for scientific analysis. 

With this motivation, in this paper wc develop a new 
approach to building mock surveys, that we dub the "onion 
universe" , which mimics the tomographic structure of pho- 
tometric surveys by decomposing the full 3D lightconc data 
structure into a set of concentric 2D all-sky maps around the 
observer. We thus propose this method as a new and efficient 
tool to exploit upcoming large photometric surveys. In par- 
ticular, this new approach allows as to make the following 
main contributions that we present in this paper: 

• Data Compression: for most of its cosmological applica- 
tions, our approach provides an effectively lossless method 
to compress simulated data by a factor ~ 1000 for arcminute 
resolution angular maps. This allows Terabyte-sized Simula^ 
tions containing tens of billions of particles to be analyzed 
in a regular laptop. 

• BAOs in the angular power spectrum: using the set of 
dark-matter density shells resulting from our decomposition 
of the N-body lightcone data, we can straighforwardly study 
the BAOs in the angular clustering of dark-matter. As an 
application of this, we assess the limit of applicability of 
linear theory (and Gaussianity of the error) as a function of 
redshift from an angular power spectrum analysis. 

• All-sky weak tensing maps: our method provides a straight- 
forward way to simulate all-sky maps of tracers of the large- 
scale structure in the light-cone from weighted combinations 
of 2D density maps. The main application of the "onion uni- 
verse" method presented in this paper is the development, 
for the first time, of an adequate all-sky simulated weak lens- 
ing map with fine angular resolution. Our map is validated 
through a comparison with theoretical predictions for the 
power spectrum over 3 decades in angular scales. We use 
this mock map to investigate effects of sampling bias in cur- 
rent weak lensing surveys from higher-order moments of the 
convergence field, and discuss the potential and limitations 
of using weak lensing as a cluster mass calibrator. 

Weak gravitational lensing by the large-scale structure 
of the universe probes density fluctuations in a wide dynam- 
ical range, from linear to highly non-linear scales. Current 
lensing observations (eg Bacon et al. 2001; Refregier et al 
2002; Jarvis et al. 2006; Masscy ct al 2007) can only sample 
the smaller scales (ie smaller than a degree) which are dom- 
inated by non-linear fluctuations. Since there is no accurate 
analytic description of the dark matter density field in the 
non-linear regime it is thus necessary to resort to numerical 
simulations to make accurate predictions of the weak lensing 
distortions. 

Simulations of weak gravitational lensing are based on 
implementations of the ray tracing technique on N-body 
simulations (eg sec Wambsganss et al 2000; Jain et al 2000; 
Vale & White 2003, and references therein). In this ap- 
proach light rays are propagated from the observer to the 



source by computing the distortion and magnification effects 
from multiple (from tens to a hundred) of equally spaced 
projectcd-mass lens planes. This approach has proven to be 
successful in measuring the lensing power spectrum which 
was found to be in agreement with the Born and Limber 
approximations (eg Jain et al 2000; Vale & White 2003). 
In addition, the lensing higher order moments induced by 
density fluctuations in the non-linear regime can also be es- 
timated (eg Jain et al 2000; White & Hu 2000; Vale & 
White 2003; White & Vale 2004) and they turn out to be 
consistent with perturbation theory results on the largest 
scales and analytic fits on intermediate scales (eg Gaztafiaga 
& Bernardeau 1998; Waerbeke et al 2000). In particu- 
lar, measurements of the variance and skewness of the lens- 
ing maps can be used to constrain the amplitude of matter 
fluctuations, erg, and matter density, f2m (Bernardeau, van 
Waerbeke & Mellier 1997; Jain & Seljak 1997). 

However, covering a wide enough dynamical range (from 
Mpc to Gpc scales) is prohibitive with current implementa- 
tions of ray tracing techniques and therefore simulations so 
far have focused on small patches of the sky (i.e few square 
degrees) , comparable to the areas covered by current surveys 
(such as GEMS, COSMOS HST, CFHTLS), well in the non- 
linear regime, where most of the lensing signal is expected 
to come from. Moreover, as we will show below, statistical 
measurements in small volume simulated surveys may be 
significantly affected by sampling bias (Hui & Gaztanaga 
1999) and cosmic variance errors are largely enhanced by 
non-Gaussianity (Scoccimarro et al 1999; Semboloni et al 
2007). For a review on the weak lensing formalism, simula- 
tions and observations see Bartelmann & Schneider (2001), 
Refregier (2003) and references therein. 

Previous work on simulating observational data in terms 
of lightcone surveys has concentrated on galaxy mocks (see 
e.g, Blaizot et al 2005; Kitzbichler & White 2007; Forero- 
Romero et al 2007 and references therein) , where finite sim- 
ulation volume, redshift discreteness, and cosmic variance 
were carefully addressed. 

This paper is organized as follows: Section 2 presents 
the simulations, the onion maps and its power spectrum. 
In Section 3, we introduce the new method to build the 
weak lensing maps from the onion shells and present several 
tests to validate the method. We also discuss different error 
estimates based on the convergence maps and a comparison 
with the halo-fit prediction. Section 4 is devoted to the study 
of non-Gaussianity, mass reconstruction and mass function. 
Finally, in Section 5 we summarize our main results and 
conclusions. 



2 ONION MAPS 

We have developed a set of large N-body simulations with 
Gadget-2 (Springel 2005) on the Marenostrum supercom- 
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Table 1. MICE N-body simulations used in this paper. 





^par 


particle mass 


Mpc/h 


number 


W^°Msun/h 


3072 


5123 


1510 


3072 


1024^ 


190 


3072 


2048^ 


24 



puter at BSC0. We shall name them MICE (MareNostrum 
- Institute de Ciencias del Espacio) simulations hereafter 0. 
In this paper we focus on a simulation run with 2048^ dark- 
matter particles in a box-size of Ltox = 3072 Mpc/h, and 
assume a flat concordance LCDM model with Qrn ~ 0.25, 

= 0.75, Qb = 0.044, = 0.95, erg = 0.8 and h = 0.7. 
The resulting particle mass is M = 2.34 X 10"Maun/h and 
the softening length used is 50 Kpc/h. Thus our simulation 
has a dynamic range close to five orders of magnitude. We 
start our run at Zi — 50 displacing particles using the Zel- 
dovich dynamics. The MICE simulation has a similar num- 
ber of particles to the Millenium simulation (Springel et al 
2005) but 6^ = 216 times more volume (and corresponding 
larger particle mass) . This makes the MICE simulation more 
adequate than the Millenium for very large scale statistical 
analyses, such as the search of the baryon acoustic scale (see 
below) and the study of very long distance effects (such as 
gravitational lensing). The drawback is the limited resolu- 
tion that is required to study smaller scales and substructure 
in galaxy size halos. In terms of volume, our simulation is 
similar to the Hubble Volume (Evrard et al 2002), but it 
has 8 times better mass resolution. Given that we have a 
relatively large particle mass, we have also used some MICE 
runs with different number of particles, given in Table [ll 
to test resolution effects. We have also varied the box size, 
ranging from Li,ox = 768 Mpc/h to Li,c,x = 7680 Mpc/h, 
to explore volume effects. Results from this analysis will be 
discussed in detail elsewhere. 

We have built a lightcone output of this simulation from 
~ 200 comoving outputs which are separated by constant 
spacing in cosmic time ( 70 Myr ). The ligtcone has been 
constructed in spherical concentric shells, each one getting 
its particles from one of the comoving outputs. The distance 
from the center of the set of spherical shells, where the ob- 
server is, to the mean radial distance of each shell is given 
by the corresponding redshift of each comoving output. In 
every shell the dark matter particles are moved using their 
peculiar velocities in order to extrapolate them into their 
lightcone positions. We have allowed a spherical shell buffer 
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Projected matter density and weak lensing maps 
from the MICE simulations are publicly available at 
|http:// www. ice. cat/mice 



to take into acount the particles that cross shell boundaries 
when moved. 

The lightcone extends out to zlc — G (i-e., a comov- 
ing distance r_tc — 6 Gpc/h from the observer) by replicat- 
ing twice the parent simulation with Ltox size along each 
cartesian direction and applying periodic boundary condi- 
tions. Note that a radial comoving distance of 3 Gpc/h cor- 
responds to z ~ 1.4, where cosmic evolution plays a signif- 
icant role for an observer at z = 0, so we do not expect 
this periodic repetition to have much impact on our anal- 
ysis. We plan to explore this issue further in future work 
by comparing results of different Ltox- Preliminary results 
indicate that L^o^ = 3072 Mpc/h is large enough for most 
applications we have explored. 

We have done several tests with the simulation out- 
put to make sure that basic clustering statistics, such as 
the power spectrum, the halo mass function and the higher 
order correlations, are consistent with previous results. As 
shown below, we have measured the baryon acoustic oscilla- 
tions (BAO) imprinted in the matter distribution at differ- 
ent redshifts or onion shells and found good agreement with 
other analyses. More details will be presented elsewhere. 



2.1 Redshift shells and BAO 

We discretize the lightcone output of the simulation in con- 
centric spherical shells of width dz given by a constant spac- 
ing in cosmic time. We use a spacing of « 70 Myr to match 
the width used for the shells in the lightcone construction 
described above. This corresponds to a dz that varies from 
dz ~ 0.005 at low z to dz ~ 0.025 at z = 1.4 (which cor- 
responds to a width of ~ 16 Mpc/h to 35 Mpc/h). This 
is probably enough for most large scale applications. Each 
onion shell is stored as a number density pixel map of given 
resolution in the Healpix format (Gorski et al 1998). Here 
we use maps with Nside ~ 2048, which pixelices the sky with 
12N^i^^ « 50 million cells of size 6pix — 1.7 arcmin size. The 
resulting onion universe decomposition of the parent simu- 
lation is shown in Fig[l] 

Fig[2][5] show four spherical shells of the onion universe 
corresponding to the projected matter density distribution 
for four different redshifts (with brighter colours represent- 
ing higher number of particles per pixel in a log scale). It is 
clearly visible a characteristic ~ 100 Mpc/h cell of filamen- 
tary structure (the so-call Cosmic Web, ie Bond et al 1996). 
The cell-size shrinks to smaller angular scales and smaller 
particle number density contrast as we move to larger red- 
shifts. By z — 0.6, there are already close to a thousand of 
these 100 Mpc/h cells in this single onion shell (i.e., Fig|5]). 
It is precisely around this ~ 100 Mpc/h scale that future 
surveys will aim at measuring the baryon acoustic oscilla- 
tions scale, vbao- The relative error involved in measuring 
tbao is roughly proportional to the inverse of the square 
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Figure 1. The onion universe: a decomposition of the lightcone 
that mimics the data structure in photometric galaxy surveys. 
The simulated universe is rendered as a discrete set of projected 
matter density shperical shells in the lightcone around the ob- 
server, i.e, at the center of the concentric spheres. 2D spherical 
shells are equally spaced in comoving time and pixelized using 
the Healpix tesselation of chosen angular resolution. For clarity, 
in this figure we only show one of the hemispheres (i.e half the 
onion universe) for several of the lowest redshift shells. 



root of the number of independent t-bao cells 

1/2 

Abao 



ArsAO 



(1) 



rsAO \ V 
where V is the sampled volume, and we have assumed Gaus- 
sian errors (with negligible shot-noise) over the first two 
BAO wiggles (see also Angulo et al 2008). Thus, for the 
onion shell at 2 = 0.6 we estimate Abao — l/VlOOO ~ 3%. 
According to this rule of thumb, we can get to 0.6% relative 
error in measuring rsAO using the whole MICE simulation 
volume, as compared to 9% with the Millenium simulation. 



2.2 Compression factor 

To build the light-cone with sufficient accuracy, we have used 
200 comoving simulation outputs. Each output takes 250 
Gbytes, so the total storage required is about 49 Terabytes. 
If we match the spatial width of the onion shells (as we have 
done) to the time lag between the outputs that are used 
to build the light-cone we will have equivalent information 
for applications that do not require angular or redshift res- 
olution better than that projected onto the pixel maps. We 
have produced 200 such Healpix maps, each occupies 201 
Megabyte, which represents a total of 39 Gigabytes. Thus, 
there is total compression factor of about 1300 when using 




Figure 2. Onion shell density map at z 2; 0.036 (this corresponds 
to a comoving distance of r = 108 ± 8 Mpc/h) 




Figure 3. Onion shell density map at 
tance r = 439 ± 9 Mpc/h) 



z ~ 0.15 (comoving dis- 
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the pixel maps instead of the full comoving outputs. Nowa- 
days, it is practically impossible to manage and share 49 
Terabytes of data in a public network: it will take more than 
a year to transfer these data over a 10 Mbit/sec connection. 
On the other hand, 39 Gbytes fits into a laptop and can 
be shared in a matter of hours. There are applications for 
which we might need the full MICE output information, but 
in many cases the pixel maps will be very useful, specially 
when we compare to observations in photometric surveys, 
as we show next. 




Figure 4. Onion shell density map at z 
tancc r = 866 ± 10 Mpc/h) 



0.30 (comoving dis- 




Figure 5. Zoom over a shell at 
r = 1589 ± 12 Mpc/h). 



0.60 (comoving distance 



2.3 Angular Power spectrum 

FigE] shows the total power per log multipole interval Pi = 
l{l + 1) Ci /2n measured from the all-sky onion maps. We 
have combined thin onion maps into dz = 0.1 slices to match 
the photo-z error expected for a survey such as the photo- 
metric SDSS (Adelman-McCarthy et al 2006) or DES (An- 
nis et al 2005f|. Symbols in the figure correspond to the Pi 
estimation in maps with redshift ranges (from top to bot- 
tom): z = OA- 0.5, 0.9 - 1.0 and 1.4 - 1.5. 

As expected, the amplitude of the fluctuations decreases 
with the depth of the slice. The BAO wiggles (around I ~ 
80 — 300), which are clearly visible, also move to smaller an- 
gular scales (higher multipoles). The scatter within multi- 
pole bins (i.e., the "intra-bin" variance) estimates the sam- 
pling variance for band limited measurements (Fosalba & 
Szapudi 2004), although this method requires appropriate 
calibration with simulations. In Section 3 we compare a sim- 
ple implementation of the intra-bin variance against other 
more standard error estimators. Shot-noise (shown as a dot- 
ted line) does not affect the Ci estimation within the range 
of scales shown. 

Linear theory predictions are shown by the continu- 
ous lines. Simulations match well the linear prediction up 
to scales of around the first BAO wiggle. On smaller scales 
(larger multipoles) non-linear effects become important. This 
effect is larger at low redshift, as expected. For the smallest 
redshift bin shown {z = 0.4 — 0.5), there is a fiattening of 
the power on scales / > 2000 indicating the virialization of 
structures. 

Fig[7]highlights the BAOs in the Ci's by normalizing to 
the linear theory prediction without baryons (i.e., same cos- 
mological parameters but fi;, = 0). The first prominent step 
at Ih — 10 corresponds to the horizon scale at the matter- 
radiation transition. The first BAO is found at around li — 
80 for z = 0.4-0.5 and h = 160 for z = 0.9-1.0. These mul- 
tipoles correspond to the projected scale for the first peak 
in the 3D power spectrum, ki ~ l\/r[z) ~ 0.07, where r(z) 



^ We take the mean depth of the survey to be 2 ~ 0.7, but this has 
a very small effect over the maps because the selection function 
is normalized to be unity in each slice 
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Figure 6. Angular power spectrum estimated from combinations 
of onion maps (symbols) of width dz = 0.10 and mean redshifts 
of 2 = 0.45 (top), z = 0.95 (middle) and z = 1.35 (bottom). 
The continuous line shows the linear theory prediction for each 
redshift. The dotted line shows the shot-noisc contribution in each 
case. 



Figure 7. Angular power spectrum normalized to a linear model 
without baryons. The continuous line corresponds to the linear 
theory prediction (with baryons) while the dashed lines marks 
the 1-sigma Gaussian errorbars. The dots correspond to mea- 
surements in the simulated maps. The long-dashed (red) line just 
shows a smooth version of the dots. 



is the comoving distance to the onion slice and the approxi- 
mation is vahd for smaU angular scales (i.e., we employ the 
Limber approximation). 

On scales where hnear theory is vahd and for i > 10 it 
is usually assumed that the Ci amplitudes follow a Gaussian 
distribution with sampling errors given by the Gaussian pre- 
diction o'q = 2Cf /{I + 1). We can test this hypothesis with 
our maps by calculating the Gaussian x^: 

E ^ ' (2) 

1=10 

where are the measured values and Cf are the linear 
theory predictions. We find = 281 and ~ 100 for 
z — 0.95 and z = 0.45 respectively. We use Imax = 230 and 
Imax = 120 to include the whole first BAO wiggle. These 
values are larger than expected for a sum of Gaussian vari- 
ables for the z — 0.95 case, and are fine for z — 0.45. The 
probabilities of this to happen are P[X^] = 70% for z = 0.45 
and P[x^] = 0.26% for z = 0.95. The later probability in- 
creases to f'[x'^] = 12% for z — 0.95 if we use Imax = 177 
(half the way through the first BAO wiggle). This indicates 
that even on these very large scales and early times there 
is correlation between different modes and one should be 
careful when doing precision forecasts using Gaussian er- 
rors. The correlation between bins is weaker if we bin the 
data in multipoles (see below). However, the problem seems 
more critical at high redshift where the BAO wiggle is bet- 



ter sampled and there are larger projection effects. A more 
detailed analysis of this effect requires more realizations and 
will be presented elsewhere. 

3 CONVERGENCE MAPS 

A basic quantity in weak lensing is the convergence field. In 
the so-called Born approximation, one integrates the lens- 
ing distortion over the unperturbed photon paths which is 
a good approximation for most applications (eg see Cooray 
& Hu 2000; Bartelmann & Schneider 2001; Vale & White 
2003; Refregier 2003). In this approximation, the conver- 
gence is just a weighted projected surface density: 

2c^ J rs a 

where 5 is the 3D matter density at radial distance r and 
angular position 6 (which is here a 2D vector) and is the 
radial position of the lensing sources|f| Out of k we can get 
all other quantities of interest such as magnification, demag- 
nification, shear, projected potential or defiection angle (see 

* Without loss of generality we will assume on writing equations 
that we live in a flat universe and that all lensing sources are at 
a fixed redshift Zs, or radial distance rg. It is straightforward to 
generalize this to an arbitrary distribution of sources, eg (Bartel- 
mann & Schneider 2001; Refregier 2003). 
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Bartelmann & Schneider 2001; Refregier 2003 for a review). 
Some of these relations are not local (they involve derivatives 
or integrals) but for an all-sky map, all these quantities are 
trivially related. In harmonic (or 2D Fourier) space, we can 
transform a convergence map into a gravitational potential 
or shear map by just multiplying the harmonic (or Fourier) 
K amplitudes by the appropriate combination of multipole 
or wave numbers. However, in the observable universe the 
boundaries of a real survey complicate these transforma- 
tions. 

We will build our convergence map by just adding the 
onion slices from the simulation with the appropriate lensing 
weight. This can be done as follows (see also Gaztanaga & 
Bernardeau 1998): 

j " ■' 

where i indicates a pixel position in the sky and j a radial 
bin (at distance rj of width drj) into which we have sliced 
the simulation as described in the previous section. If we 
indicate by Nij the number of particles in pixel i from onion 
slice j, we have: 

<5(i,j) = ^-l (5) 
where p =< p{i,j) > and 



where Af2 is the area of each pixel. 

FigH] shows images of the resulting maps for lensing 
sources at = 1- Note how despite the large volume pro- 
jected there is still considerable structure in these maps. In 
particular, on scales of a few square degrees there is a large 
variation from place to place in the maps. This indicates 
that sampling variance is important. Current weak lensing 
surveys, such as COSMOS (Massey et al 2007), or lensing 
simulations, only expand scales of the order of a few square 
degrees. Our convergence maps show that sampling variance 
is quite large on such small scales and it is unlikely that cur- 
rent data could represent a fair sample of the universe. In 
section 4 and in the Conclusions we will show some more 
quantitative consequences of this note. 

3.1 Validation of the map 

We can validate the convergence map by comparing its power 
spectrum and higher orders to theoretical expectations. First 
we will focus on a comparison with linear theory and a con- 
sistency test for the angular power spectrum. 

Based on its definition in Eq|3] we would expect the 
angular power spectrum of the convergence to yield: 

C.{n)^'J^ j drPik,.)^Il^ (7) 

where P{k, z) is the 3D density power spectrum in the sim- 
ulation at redshift z (corresponding to the radial coordinate 
r = r(«) in the integral) evaluated at fc = //r in the small 
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angle (Limber) approximation, valid for I > 10 within a few 
percent accuracy (see e.g, Vale & White 2003). In terms of 
discrete onion shells, this translates into: 

= E Pil/r.,^.)^-^'^ (8) 

Fig|9] shows a comparison of the angular power spec- 
trum in the above prediction (continuous line through the 
symbols) to the power spectrum directly measured from 
the convergence maps (symbols with errorbars). The power 
spectrum has been binned in adjacent multipoles with bin- 
width Al = 20. The errorbars indicate the scatter of power 
within a bin. For the prediction in Eq|8] we have used the 
actual 3D power spectrum measured from all particles in 
the comoving outputs of the MICE simulation at the corre- 
sponding redshifts. The dot-dashed line uses the density in 
a PM grid of 2048^ to estimate the same P{k,z). In both 
cases, this is just an approximation because P{k, z) should 
be the power spectrum in the lightcone. But the difference 
is small because the redshift shell is quite narrow. It is also 
an approximation in the sense that P{k, z) in the whole box 
could be different to the power in a particular onion shell 
(redshift bin), due to sampling variance. 

The measured spectrum agrees with the linear predic- 
tion on the largest scales (I < 100), as expected. The agree- 
ment in shape and amplitude with the prediction validates 
the convergence maps on the largest scales. On all scales the 
agreement with Eq[8]is excellent, indicating that the way we 
have built the convergence maps is a good approximation to 
the true map. It also indicates that statistically, inhomo- 
geneities in the radial direction do not affect the projection, 
which was assumed to get to Eq|8l 

The deviations at I > 1000 when using the PM grid 
P{k, z) (dot-dashed line) show that the power spectrum on 
those scales come from structure on scales smaller than the 
cell size used for the PM grid ~ 1.5 Mpc/h. Thus, one needs 
higher resolution PM simulations (or treePM algorithms for 
a given PM grid size) to model the larger multipoles. 

In Fig llOl we show the cumulative contribution in Eq[8] 
from 2 — Zmin to z = 2s . For multipoles I > 100 the relative 
contributions are quite flat with 2%, 20%, 50% and 85% 
contribution coming for sources at Zmin > 0.8, Zmin > 0.6, 
Zmin > 0.4 and Zmin > 0.2, respectively. All redshifts seem 
to add power in the non-linear multipoles / > 1000. Note 
how the lowest multipoles have a contribution > 50% from 
z < 0.2, due to local structures. 

3.2 Error comparison 

The convergence map with sources at 2s = 1 shown in 
FigO nominally covers a volume, Vn ~ 4/37rr(2s = 1)"' — 
58 Gpc^/h^ with r{zs = 1) ~ 2400 Mpc/h, which is twice 
the volume of the parent MICE simulation, Vp = Lf^^ = 
(3072 Mpc/h)^ ~ 29 Gpc^/h^. This is possible because we 
have replicated twice the box in each cartesian axis to get 




Multipole 

Figure 9. Angular power spectrum in the convergence maps 
(symbols with errorbars) as compared to linear theory (dashed 
line) and predictions from the full measured 3D power spectrum 
in Eq[7](line that goes through the symbols). The dot-dashed line 
uses the same prediction with the 3D power spectrum measured 
only using a 2048'' particle-mesh (PM) grid. 




10 100 103 

Multipole 



Figure 10. The top panel shows as dot-dashed, long-dashed, 
short-dashed and dotted lines the cumulative contributions to 
the sum in Ec[|8] starting from z — ^min 

= 0.8,0.6,0.5 and 0.2 
respectively (as labeled in the figure) and integrated to 2 = Zs = 
1. The bottom panel shows the ratio of these quantities to the 
total contribution. 
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to z ~ 6 (see also Section 2). However, note that to reach to 
Zs — 1 each sky octant (1/8 of the full sky) is truly indepen- 
dent in the way we have built the lightcone. Firstly because 
each sky octant uses the same replica but with the observer 
placed in a different box position, so the same structures 
are seen differently (i.e., they are seen from a different angle 
and distance). Secondly, because of the shape of the lensing 
efficiency, the shells closer to Zs = 1 (where the shell- volume 
is larger) give a negligible weak lensing signal. We can de- 
fine an effective volume 14 sampled by the convergence map 
as the volume weighted by the lensing efficiency function 
(renormalized to be one at the peak efficiency). We find 
that in fact 14 ~ V„/2 ~ 29 Gpc^/h^. We thus conclude 
that the convergence map to Zs — 1 samples well the full 
parent volume without significant repetition. 

We can therefore split the sky in fO equal area disjoint 
regions (a partition each f 0% of the sky) and measure Ci in 
each of the fO regions using a fast 2-point estimator, SpICE 
(Szapudi et al 200fa; Szapudi et al 200fb). We use the vari- 
ance in the Ci from each region to have a direct estimate of 
the errors for a survey that covers a fraction fsky = O.f of 
the sky. We call this the "sub-sample" (SUB) error. To avoid 
the covariance between bins, we follow (Cabre et al 2007) 
and bin adjacent Ci estimations in bins of with Al = 20 
for / < too and Al = 40 for I > 100. These binwidhts ren- 
der the covaraince matrix effectively diagonal for the binned 
Ci's (see Cabre et al 2007 for details). We can also use a 
variance estimator based on the rms dispersion of individual 
Ci amplitudes within in a given bin to estimate the error in 
that bin. We shall call this "intra-bin" (IB) error estimation 
(Fosalba & Szapudi 2004). For Gaussian distributed ampli- 
tudes this is a good estimate, because there is no correlation 
between adjacent multipoles. 

We will compare the above errors with the traditional 
Gaussian estimate of the sampling variance (SV): 

where Al — Ntin is the number of multipoles in each bin. 

Fig lffI compares the different estimates for the relative 
errors in Ci{k,). On scales I > fOOO there is a good general 
agreement. SV errors seem to underestimate SUB errors by 
about ~ 50% between I = 500 and I = 2000, but they yield 
compatible results otherwise . The IB estimator does well for 
I > 500 but can be a up to factor of 4 too large for I < 500. 

3.3 Power spectrum: mass resolution and 
shot-noise 

Fig lf 21 shows the power spectrum in the convergence maps 
with different mass resolution (but same Li,ox, see Table [T]). 
There are two different contributions to the effects shown 
here. One is the possible difference due to mass resolution 
and the other is due to the finite particle density, which 
results in a different shot-noise contribution. 




Multipof 



Figure 11. Comparison of relative errors in C;(k) for 10% of 
the sky: a) variance from 10 subsamples in all-sky convergence 
map (continuous line, SUB) b) sampling variance from Gaussian 
statistics (dashed line, SV) and c) intra-bin variance (long-dashed, 
IB). 
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Figure 12. Raw power spectrum estimated from convergence 
maps in simulations (symbols) with Zs = 1 and for different reso- 
lutions (shown in Table [TJ. Continuous lines show the linear and 
non-linear halo-fit predictions which are obtained by replacing 
P(k,z) in Eq[7]by the corresponding 3D power spectrum. 

To correct for shot noise we subtract from the mea- 
sured Ci a term given by Eq[8]with P{k,z) = t/N, where 
N is the mean galaxy density at each redshift. The corrected 
spectrum is shown in Fig[T3l The results for the two higher 
resolutions agree well up to Z ~ f 000, and roughly follow the 
halo-fit prescription. The power spectrum estimated from 
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Figure 13. Same as Fig |12l but corrected for shot-noise. The 
correction in each case is shown by the dotted lines. 

the low resolution map, containing 512'^ particles (N512) 
and mass resolution M ~ 1.5 x lO^'^Msun/h), closely fol- 
lows the linear (rather than the non-linear) prediction. This 
could be due to the lack of small halos and it is also appar- 
ent when we compare measurements of the 3D power spec- 
trum in outputs from the low resolution measured in the 
low resolution 512^ MICE comoving simulations outputs to 
the corresponding linear P(k, z) predictions. The inability of 
the low resolution simulation to reproduce non-linear effects 
reflect the well known fact (e.g., see halofit in Smith et al 
2003) that the power on non-linear scales is dominated by 
the internal structure of halos with mass smaller than a few 
times lO^^Maun/h, i.e., comparable to the particle mass of 
the low resolution simulation. 

3.4 Power spectrum: halo-fit 

Fig |14l shows the relative difference between the measure- 
ments in the (all-sky) simulations and the halo-fit predic- 
tion. This prediction is obtained by using the halo-fit model 
P{k, z) in Eq[71 The dashed lines show the dispersion in 10 
subsamples (each 10% of the sky). The halo-fit model only 
seems to work within 5% accuracy up to multipoles I < 1000. 
Deviations at smaller scales are significant (up to 30%). This 
is the case even for a survey which is only 10% of the sky 
(i.e., fsky = 0.1), shown as dashed lines in Fig |14l Note that 
contrary to Fig llll we have not binned the data here. This 
explains the large discrepancy between the Gaussian errors 
(dotted lines) and the subsample errors (short-dashed lines). 
Correlation between bins is strong resulting in smaller di- 
agonal errorbars, but larger covariance for the subsample 
errors. Surprisingly, the Gaussian error prediction for the 
all-sky simulations (long dashed lines) is similar to the mea- 
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Figure 14. Points show the angular power spectrum of the all- 
sky convergence maps relative to the non-linear halo fit model. 
The continuous line is a smoothed version of the dots. The short- 
dashed line indicates the errors in 10 subsamples with 10% of 
the sky. The dotted (long-dashed) lines show the corresponding 
Gaussian error for 10% (100%) of the sky. 

sured diagonal error in sub-samples of fsky =0.1 (see also 
Cabre et al 2007 for a related discussion in cross-correlation 
analyses) . 

4 NON-GAUSSIANITY AND PROJECTION 
EFFECTS 

4.1 Moments 

It has now been well established that the p-order cumulants 
of the projected local density field {k^)^ are expected to 
behave as 

{.'')^ = S,{.y (10) 

with p = 3, 4..., on large scales (see Beranardeau et al 2002 
and references therein). The Sp parameters, quantify the 
departure from Gaussian behaviour, and the variance (k^)^ 
can be obtained from the power spectrum above. The mea- 
surement of the gravitational weak shear induced by the 
large scale structures in deep galaxy catalogs will reveal this 
correlation properties of the projected mass, at the level 
of the two-point function (Blandford et al. 1991, Miralda- 
Escude, 1991, Villumsen 1996, Jain & Seljak 1996, Kaiser 
1995) or for higher orders (Bernardeau, van Wearbeke & 
Mellier 1996, Gaztafiaga & Bernardeau 1998). 

Figure lTSl shows the above moments for the convergence 
maps, smoothed with a Gaussian window of size 9. Results 
are in good agreement with the analytic predictions men- 
tioned above. More details of this comparison will be given 
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Figure 15. Squares show the variance (bottom), skewness (mid- 
dle) and kurtosis (top) measured in the all-sky convergence map 
with Zg = 1. Triangles show that the mean in 10* small subsam- 
ples of 1.6 sqr. degrees that match the size of the COSMOS HST 
field is biased with respect the result using all the sky. The larger 
errorbars represent the 1-sigma variance in the COSMOS HST 
subsamples. The smaller errorbars correspond to fsky = 0.1. 
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Figure 16. Squares show the mean ratio of reconstructed versus 
true density fluctuations in the convergence maps. Continuous 
lines represents the median of the distribution (50% percentile). 
The short-dashed (long-dashed) line show 25% (16%) and 75% 
(84%) percentiles. Each panel corresponds to a different redshift, 
as labeled in the plots. 



elsewhere. Here we just want to stress that the maps are 
strongly non-Gaussian and we want to focus on the error 
estimation. 

The smaller errorbars in Figure[T5]correspond to fsky = 
0.1, while the larger errorbars correspond to a field size of 
1.6 square degrees, comparable to the HST COSMOS field 
(Massey et al 2007). Squares correspond to the mean in the 
all-sky map, which agrees very well with the mean of the 10 
fsky =0.1 subsamples. This is not the case for the mean 
in the (10*) COSMOS-sized subsamples, which is severely 
biased. We find a 10% bias for the variance and skewness 
at 2', and this relative bias increases by a factor of 2 for 
10' and a factor of 4 for 30' scales. This sampling bias is 
common when one has large fluctuations at the scale of the 
survey (eg Hui & Gaztanaga 1999), as is the case here (e.g., 
see also FigOfor a visual impression of this effect). 

4.2 Mass reconstruction 

We could in principle use the weak lensing signal to map 
the 3D mass distribution in the universe (e.g., Massey et al 
2007) or we can use it for the more modest task of calibrating 
the mass of known clusters with the idea of estimating the 
cluster mass function (see e.g.. White, Waerbeke & Mackey 
2002). Here we want to use the large statistics in our huge 
volume simulation to study how important projection effects 
could be for this mass calibration. 



At each sky pixel i we flnd the redshift onion shell j 
where the contribution to K(i) in Eq|4] is maximum. This 
produces an array of 3D pixels which are potential sites 
for large overdensities (e.g., those one can typically associate 
to clusters in galaxy surveys) which we want to use for mass 
reconstruction. We will assume that both j and Za are known 
and we will use the total convergence K{i) to estimate the 
overdensity at as: 

This procedure mimics a simple linear mass reconstruction 
method which assumes that the measured convergence K{i) 
is dominated by the overdensity of pixel (i, j) (the maximum 
along the line-of-sight). 

Square symbols in Fig llGl show the mean ratio between 
the reconstructed 5{k), given by Ea llll and the true mean 
fluctuation 5{i,j) over all pixel positions i in the simulation, 
as a function of the true fluctuation. Each panel corresponds 
to a different redshift (i.e., a different value of j in Eq llip . 
The long-dashed lines correspond to the 1-a scatter in the 
reconstruction. As can be seen in the flgure, the mean re- 
construction is basically unbiased but there is quite a large 
scatter which increases as we decrease the size of the true 
fluctuation. 

The top panel in Fig |17l shows cumulative histograms 
of the above reconstruction, where we have converted the 
fluctuations 5 into mass using M — {1 + 5)pdV, where p is 
the mean pixel density and dV is the pixel volume. Note that 



© 0000 RAS, MNRAS 000, 000-000 



12 



Fosalha et al. 



10-' 

102 



z = 0.5 



> 2x10'-' 




0.6 0.8 1 1. 

M/M„.„„. 



1.4 



1.6 

S 

1.4 
\ 1.2 

0.8 
J 0.6 



















^ E ffi ^ 








m 












m 












m 






o : 






m 










S 




b ~ 






z = 0.b^O,l 


K4= 1233«: 


1 - 
OJ — 






, , , , 1 


1 





10 100 
Mass (xlO'2 solar) 



Figure 17. Top: Histogram showing tlie ratio M/Mf^appat of 
true versus recovered convergence mass, for pixels in the simu- 
lation with true mass M above xlO^^Msun/h (outer histogram), 
5 X lO^^Msun/h (middle histogram) and lO^'^Msun/h (inner his- 
togram). The mean recovered mass is not biased, but the distri- 
bution is quite broad and non-Gaussian. Bottom: implication of 
the above histograms for the recovered mass function. The ratio 
of the recovered over the true cumulative number of pixels above 
a given mass is shown as a function of the true mass. Errorbars 
represent 2-a rms dispersion when we split the simulation in 10 
pieces of fsky =0.1 each. 



the distribution is not Gaussian. These results agree, at least 
in a qualitative way, with more detailed studies over smaller 
simulations (e.g.. White, Waerbeke & Mackey 2002). 



4.3 Implications for the mass function 

The large scatter in the above mass reconstruction has im- 
portant implications for estimating the (cluster) mass func- 
tion. This is illustrated by the bottom panel of Fig |17l which 
shows the ratio of the recovered versus the true cumulative 
mass function ai z — 0.5 ± 0.1. We chose this redshift be- 
cause it is effectively where the convergence window function 
is maximum for sources a,t Zs — 1. We also show 2-a error- 
bars from the variance in 10 subsamples with fsky ~ 0.1. As 
can be seen in the figure, there is a significant bias, given 
the errorbars, in the mass function. Deviations can be as 
large as ~ 50%, with an excess density at the high mass 
end (M > 10^"^ Msun/h) and a smaller deficit at the lower 
mass end. This is expected because there is a larger num- 
ber density of smaller mass overdensities, which results in a 
greater excess of smaller mass objects that scatter into the 
large mass bins. 



5 SUMMARY & CONCLUSIONS 

Radial shells are a natural and convenient decomposition 
of the data volume to exploit large astronomical surveys. 
Both because of the hmitations in our ability to measure 
precise redshifts for all objects in galaxy surveys (especially 
at z > 1), and also because it is a natural way to split the 
survey data to study evolution or avoid redshift space distor- 
tions. Here we have presented a new generation of very large 
scale N-body simulations developed at the Marenostrum su- 
percomputer using the GADGET-2 code. We have named 
them the Marenostrum Institut de Ciencies de I'Espai or 
MICE simulations. In this first paper we have focused on a 
simulation that contains almost 10^*^ particles in a cubical 
box of 3 Gpc/h on a side, what delivers good enough mass 
resolution (2.34 x 10^^ Msun/h) for studying clustering statis- 
tics of the large scale structure and cover a dynamical range 
of five orders of magnitude: from Gpcs to tens of Kpcs. This 
allows us to sample from the largest (linear) scales to very 
small (non-linear) structures. 

Using this MICE simulation we have built an all-sky 
lightcone output that extends to high redshift by replicat- 
ing the parent simulation box Lbox- We have shown that, 
thanks to the way we build the lightcone, the effect of re- 
peated structures at distances r > Ltox from the observer 
is effectively negligible in clustering measurements. In order 
to mimic the onion-like structure of real data from galaxy 
surveys we have compressed the lightcone data into a set 
of radial shell maps of given redshift resolution. These all- 
sky angular maps have then been pixelized using the conve- 
nient Healpix tesselation with high spatial resolution. Our 
approach provides an effectively lossless method to compress 
simulated data by a factor ~ 1000 for arcminute resolution 
angular maps. This allows Terabyte-sized simulations con- 
taining tens of billions of particles to be analyzed in a regular 
laptop. 

We have presented two main applications of the onion 
maps for the study of large-scale sctructure in the universe: 

• the study of BAOs in angular maps of the dark-matter 
distribution in the lightcone 

• the construction, for the first time, of an adequate all- 
sky simulated weak lensing map with fine angular resolution. 

The onion maps we have generated are large enough 
to detect the BAO scale with a precision better than 1%. 
We have presented the angular power spectrum of the maps 
(see Figs l6l7|l and discussed how accurate is the Gaussian 
estimation of sampling errors in the presence of non-linear 
effects. 

In Section 3, we have used the onion shells to build a 
new set of all-sky lensing maps. These maps are validated by 
comparing its power spectrum and higher orders in the map 
to theoretical expectations. In Fig llOl we show the relative 
contribution of each redshift shell to the total convergence 
power. Because we simulate the entire celestial sphere we can 
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compare theoretical prescriptions for error estimation with 
an error based on the rms dispersion over sky patches (sub- 
samples). As summarized in Fig llll we find that Gaussian 
errors tend to underestimate the true errors in the power 
spectrum by about ~ 50% at non-linear scales in multipole 
bins between I = 500 and / = 2000, even if we use broad 
bins to minimize covariance between adjacent multipoles. A 
comparison of the convergence power spectrum to the ana- 
lytic halofit predictions yields discrepancies of order 30% on 
highly non-linear scales, I > 1000 (see Fig |12lli|) . 

We find that current lensing surveys (such as Massey 
et al 2007) might yield biased estimates of the clustering 
statistics since measurements are subject to large sampling 
variance. In other words, these small surveys do not repre- 
sent a fair sample of the universe. This has been quantified 
in the variance and higher-orders of the convergence field 
(see Fig|15p by randomly sampling the all-sky lensing map 
with ~ 10** COSMOS-sized surveys. For the variance, we 
estimate a ~ 10% bias and ~ 40% errors at 2', and the bias 
increases to 25% and 50% at 10' and 30' scales. We also find 
comparable relative biases for the skewness as a function of 
scale. 

Our estimate of the error is significantly larger than the 
value reported in Massey et al (2007) for the COSMOS sur- 
vey. The origin of this discrepancy might be the fact that 
Massey et al (2007) compute the variance using subsamples 
of the COSMOS data, and thus they do not include sample 
variance at the scale of the survey. Consequently, cosmo- 
logical constraints on asil'^''^* based on this estimate of the 
variance are expected to be also biased low significantly and 
the error could be underestimated by as much as as a factor 
of 2. 

We have also measured the degree of non-Gaussianity 
in lensing maps induced by non- linear growth. Higher or- 
der moments in the convergence field, shown in Fig |15l are 
compatible with theoretical expectations. In particular, they 
match well the amplitudes of the hierarchical scaling ex- 
pected from non-linear gravitational clustering (Beranardeau 
et al 2002). In the case of weak lensing, these amplitudes are 
also strongly dependent on cosmological parameters Bernardeau, 
van Waerbeke & Mellier 1997; Gaztanaga & Bernardeau 
1998. 

Finally, we have presented a mass calibration proce- 
dure using lensing maps. In Fig |17l we illustrate how well we 
expect to recover mass estimates based on all-sky conver- 
gence maps. Upcoming wide surveys, such as DES (Annis 
et al 2005), plan to calibrate the cluster mass function us- 
ing the weak lensing information. We have shown here that 
this approach is a promising tool for calibrating masses, but 
it needs to be corrected from systematic biases that arise 
because of the large scatter induced by projection effects. 
Further work needs to be done to check the impact of these 
considerations in specific mass callibrationmethods, such as 
the one recently presented in Johnston et al. 2007. 



Note added in proof: after our paper was submitted, 
other papers appeared presenting applications of N-body 
simulations for CMB lensing analyses (Das & Bode 2008; 
Carbone et al. 2008). We note that the method introduced 
in our paper can be straightforwardly applied to CMB lens- 
ing. We will present this application elsehwere. 
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